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1. Introduction 

It is a characteristic feature of non-Hermitian but PT-symmetric Hamiltonians that, 
as an external parameter is varied, pairs of real eigenvalues coalesce at an exceptional 
point (EP) and turn into complex conjugate pairs (Bender and Boettcher 1 , Moiseyev 
[2]). Recent investigations of non-linear PT-symmetric Hamiltonians modelling Bose- 
Einstein condensates in double-well potentials with sink and source terms of atoms 
at the respective wells have revealed a specific spectral behaviour which implies 
apparently new spectral singularities [3HS]- Pairs of real eigenvalues merge at an 
EP when the strength of the non-Hermiticity is increased, but no pairs of complex 
conjugate eigenvalues appear beyond the EP. Rather, such pairs are born as if "out of 
nowhere" before the EP, that is for smaller values of the loss/gain term, where they 
bifurcate from the branch of the ground state eigenvalue. 

The reason for this unusual behaviour is the non-analyticity of the non-linear term 
ij) in the underlying Gross-Pitaevskii equation. Numerical calculations for Bose- 
Einstein condensates in PT-symmetric double- well potentials have shown [7HS] that by 
appropriate analytic continuation of the non-linear term pairs of complex conjugate 
eigenstates indeed arise at the branch point. Moreover, in an analytically solvable 
toy model for a PT-symmetric two-mode Bose-Einstein condensate Graefe [B] used 
a particular analytic extension of eigenvalues and eigenstates. This way it has been 
demonstrated that the "prematurely" born pairs of complex conjugate eigenvalues 
emerge from an additional pair of real eigenvalues. Therefore three real eigenvalues 
effectively coalesce at what we shall term a triple point. 

The features of a non-linear term of the form ip have become very relevant for 
the study of 'PT-symmetric physical systems since this type of non-linearity not only 
appears in the mean-field description of Bose-Einstein condensates. PT-symmetric 
optical setups of coupled dual waveguides including a Kerr non-linearity show exactly 
the same structure in the underlying equations. They have been used to demonstrate 
the existence of uni-directional structures [10] or the presence of solitons in loss/gain 
media [TTUT6] . The remarkable success of realising PT symmetry and PT-symmetry 
breaking experimentally emphasises the relevance of these optical setups [17l[T8] . The 
additional properties of the combination with a non-linearity as the possibility of uni- 
directionality 10 might lead to new technical devices. In quantum mechanics non- 
linear PT-symmetric systems have been discussed in model potentials by Musslimani 
et al. [19] . and for Bose-Einstein condensates described in a two- mode approximation 
[SHE] or by solving the mean-field limit of the Gross-Pitaevskii equation in position 
space [7HH]. 

Potentials consisting only of delta-functions often provide a very simple access 
to solutions, in many cases they can be obtained analytically. For example, the 
spontaneous symmetry breaking in double-well structures has been investigated 
analytically in the limit of infinitely narrow potential wells described by delta-functions 
by Mayteevarunyoo et al. [20] • Rapedius and Korsch investigated the decay of Bose- 
Einstein condensates in a double-delta-setup [H] . Bifurcations of stationary solutions 
have been studied by Witthaut et al. [22] for Bose-Einstein condensates in a delta- 
comb. With the exception of the study of Bose-Einstein condensates in a PT- 
symmetric double-delta-trap [7l[8] , the discussion of PT-symmetry in delta-potentials 
has been restricted to linear quantum mechanics. Jakubsky and Znojil [23] considered 
analytical solutions of a particle in an infinitely high square well containing two purely 
imaginary delta- functions in a PT-symmetric arrangement. The influence of a pair of 
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Figure 1. Real parts (a) and imaginary parts (b) of the eigenvalues re of the 
full non-linear equation ifTJ as functions of the loss/gain parameter 7 for a = 1.2 
for the three values of the non-linearity parameter (in (a) from bottom to top) 
g = 0.0, 0.15 and 1.0. Solid lines denote purely real eigenvalues, dashed lines 
complex conjugate eigenvalues. For g 7^ the complex conjugate eigenvalues 
bifurcate from the ground state branch before the branch point where the two 
real solutions coalesce. 



■pT-symmetric delta-functions on the bound and scattering states of a single real delta- 
function has been investigated by Jones [23]. The spectral properties, in particular 
bound states and spectral singularities, in non-Hermitian delta-potentials have been 
studied by Mostafazadeh et al. [25rl27] . 

It is the purpose of the present paper to study the new types of spectral 
singularities in more detail. We do this in Section [2] by solving numerically an 
analytically continued version of the Gross-Pitaevskii equation for the model of a 
Bose-Einstein condensate in a "PT-symmetric double delta potential. We find that 
the triple point can exhibit the behaviour of a second-order (EP2) or a third-order 
(EP3) exceptional point depending on the particular encirclement. In a matrix model 
presented in Section 3 this behaviour is described in analytic terms. 

2. Complexification of the Gross-Pitaevskii equation 

2.1. Model with delta-function traps 

The model of a Bose-Einstein condensate trapped in two PT-symmetric potential 
wells represented by delta- functions has been introduced in [7J|S]. Here we only 
recall the essential points which are relevant to the present discussion. In the 
mean-field description the non-linear Gross-Pitaevskii equation to be solved reads 
in dimensionless units 

- ip"{x) - [(1 + vy)S(x + a/2) + (1 - vy)S(x - a/2)} ^(x) 

-<#(Z)|V(*) = ~K^{x). (1) 

Here g is the non-linear interaction strength, the real- valued parameter 7 determines 
the strength of the gain and loss terms in the two delta wells located at ±a/2. 
(The relation to physical quantities can be found in [?1[5] as well as details on the 
numerics.) Since we are interested in the complete eigenvalue spectrum, including 
complex eigenvalues, solutions are searched with n G C. Figure [1] recapitulates the 
results for the eigenvalues k for different values of the non-linearity parameter g. It can 
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Figure 2. Real parts (a) and imaginary parts (b) of the eigenvalues n for the 
states appearing at g = 1 and g = 0.15. Using the substitution |i/>(x)| 2 — > 
ip(x)ip(—x) the analytic continuation beyond the branch point 7b p of the two real 
"PT-symmetric eigenvalue states is obtained. 



be seen that for every value of g a pair of real eigenvalues exists up to some value 7b p , 
where they coalesce at a branch point. Pairs of complex conjugate eigenvalues emerge 
at critical values 7 cr < 7b p , where 7 cr decreases for increasing g. For sufficiently large 
non-linearity these solutions already appear for 7 = 0. Thus one has ranges of 7 where 
two real and two complex conjugate eigenvalues coexist. The wave functions show 
that those belonging to the real eigenvalues are themselves PT-symmetric, i.e. they 
are eigenfunctions of VT, whereas the two states with complex eigenvalues k are not. 

Figure [1] reveals that the branching-off of the real eigenvalue changes continuously 
from the branch for g = 0. The wave functions, however, show a non-uniform limit as 
g tends to zero: for every value g > their asymptotic form is given by 1/ cosh(/ca;) , 
while at g = it is exp(±Ka;). The non- uniform behaviour of g — > will be addressed 
in a few cases below. 

Pairs of complex conjugate eigenvalues emerge at the branch points only when 
the non-analytic term |^(a:)| 2 in (fl]) is continued beyond the branch points. This is 
illustrated in figure [2] For the analytic continuation we exploit the 'PT-symmetry 
of the wave functions corresponding to the real eigenvalues, i.e. ip*( x ) = VK — x )- 
Therefore, when approaching the branch point we replace the non-linear term g\ip(x) | 2 
for the PT-symmetric states by gijj(x)ip(—x). This function can be continued 
analytically. We note that in the numerical calculation the additional condition 
J ip(x)ip(—x)alx = 1 must be enforced to fix the phase of the non-linearity in the 
VT broken regime. 

The substitution of |^(a;)| 2 by ip{x)tp{— x) allowing analytic continuation of the 
square modulus of PT-symmetric eigenstates shows that the number of solutions does 
not change at the branch point 7b p . Yet, at 7 < 7 cr < 7b p the number of states 
seem to be less. In fact, the pair of complex eigenvalues bifurcating from the ground 
state at 7 cr seem to come about without real-valued precursors. However, employing a 
specific analytic extension of the Gross-Pitaevskii equation these precursors and their 
corresponding eigenstates can indeed be found. 
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Figure 3. Real parts Re(ft) (a) and imaginary parts Im(ft) (b) of the spectrum 
after the analytical extension for g = 1. The red solid and green dashed 
lines represent the "PT-symmetric and "PT-broken states already known without 
any analytic continuation. The new states only visible with the full analytic 
continuation arc marked by blue dotted lines. With these continued states the 
number of states does no longer change with 7. 



2.2. Analytic continuation to 7 < 7 cr of the states breaking the VT -symmetry 

The procedure [28] is to express all complex quantities by two real-valued functions, 
either amplitude and phase, or real and imaginary parts. The non-analytic 
Gross-Pitaevskii equation then decomposes into two coupled real-valued differential 
equations which can be continued analytically. 

We write for the wave function, the double well potential and the eigenvalue 
ip = ip T + i^i, V — V T + iVi and k = K T + in- u respectively. Inserting into the Gross- 
Pitaevskii equation and sorting the real and imaginary contributions leads to 

- V" + V x ip T - Vitpi - g(tf + ipf)ip T = K T ip T - KiVi , (2a) 

+ + ViV-r - flW>r + = M'i + KiV'r ■ (26) 

The analytical extension consists in allowing the real and imaginary parts of the wave 
function and the eigenvalue to become complex again: 

ip z = ip„ + i^ri , ipi = ipu + iV'u , (3 a) 

— — fcrr ~\~ i/ijri ; = ^ir ~\~ i^ii ■ (36) 

When plotting the eigenvalues we recombine the 4 real quantities K rr , K r j, K; r , Ka into 
the real and imaginary parts of k, i. e. 

k = (K rr — Kji) + i(K ri + K ir ) = Re(K) + ilm(K) . (4) 

For the value of the non-linearity parameter g — 1 the results for Rc(k) and Im(K) 
are shown in figure [3] For 7 < 7 cr there are indeed two new branches of real- valued 
k. They merge at the triple point and give rise to the complex conjugate pair of 
eigenvalues. Moreover, these two solutions are present even at 7 = 0, where they refer 
to two degenerate eigenvalues. 

It is remarkable that the spectra of the two additional states obtained by the 
analytic continuation become very similar to the spectra of the original 'PT-symmetric 
solutions when g is decreased. This is demonstrated in figure |4l where the spectra 
shown should be compared with those in figure[T] There is, however, one difference: the 
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Figure 4. Eigenvalue spectra for g = 0.15 (red solid lines), g = 0.10 (green 
dashed lines), and g = 0.05 (blue dotted lines). For comparison the spectrum 
without the non-linear term (g = 0) is also displayed . 

two additional solutions are degenerate at 7 = for any g ^ 0, while the eigenvalues 
of the two original 'PT-symmetric states always have different values at this point. 

For g > the wave functions of the analytically continued states shown in figure 0] 
differ of course from those of the original real eigenvalues where only k„ features. It 
is of interest to have a look at these wave functions. They are illustrated for g = 1 in 
figure [5] for a few values of < 7 < 7 cr . Note that the functions are eigenfunctions 
of VT throughout in that their real parts are symmetric and their imaginary parts 
antisymmetric. Note further that there is a non-vanishing imaginary part ip T i due to 
the analytic continuation even at 7 = 0. At this point the two wave functions are 
associated with the same energy value as Ka — 0, i.e. there is a genuine degeneracy. 
(This in in contrast to the wave functions of the other two energies appearing for 
7 = (see figure [3]) ) . The imaginary part i/'ri tends to zero when 7 cr is approached 
while another imaginary part ipir ~ n °t due to the analytical continuation - emerges. 
Eventually, that is at 7 cr , the two wave functions become equal to themselves and to 
the third (original) state with which they coalesce at the triple point. From this point 
onwards the states are found without the complex extension, in fact the parts invoked 
from the extension vanish. Note also, that the additional real part ipu switches on and 
off in the interval [0,7 cr ] similar to the behaviour of na. 

2.3. Exceptional point behaviour 

As pointed out above the three eigenfunctions are identical at the triple point. 
Therefore the question arises whether or not the triple point is a third-order 
exceptional point (EP3). The way to test the exceptional point behaviour numerically 
is to extend one parameter into the complex plane, encircle the point of coalescence 
and check the interchange of the eigenvalues and states. Figure [6] shows the result for 
a non-linearity of g = 1 when the triple point at 7 cr ~ 0.308 is encircled in the complex 
extended 7 plane. One clearly sees that only the two complex eigenvalues interchange 
while the ground state remains unaffected. The same holds for the corresponding 
eigenfunctions. It thus appears that we deal with an EP2 corresponding to the merging 
of the two new eigenvalues. 

However, when we allow an arbitrarily small asymmetry for the real part of the 
double well potential by considering V asym = A [S(x — a/2) — S(x + a/2)], and encircle 



S = 0.15 

g = 0.10 - 




Spectral singularities in VT -symmetric Bose-Einstein condensates 



7 




Figure 5. Wave functions of the analytically extended solutions for a few values 
of 7 and g = 1. The right (left) column refers to the lower (upper) branch of the 
eigenvalues. 



now the branch point around A — in the complex extended A plane, figure [7] clearly 
points to an EP3. We stress that any other perturbation would produce the same 
result: the apparent EP2 which does not involve the third state becomes a true EP3 
nesting together all three states as soon as the problem becomes disturbed, no matter 
of what nature and how small the disturbance may be. This behaviour is similar to that 
of a study of a simple matrix model by Demange and Graefe [29| where three coalescing 
eigenvectors can exhibit both the behaviour of second- or third-order exceptional point 
depending on the choice of parameters. The same behaviour has also been found in 
the investigation of pitchfork bifurcations occurring in Bose-Einstein condensates with 
dipolar interaction [3Tj] . 

To summarise: the numerical solution of the 'PT-symmetric Gross-Pitaevskii 
equation (JXJ) brings about spectral properties calling for deeper insight of the physical 
system. In the following section the precise nature of these spectral singularities will 
be clarified in analytic terms by means of a matrix model simulating the properties of 
the non-linear physical system. 
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Figure 6. Exceptional point behaviour: circle with complex 7 = 0.308 + 0.04e 1 *', 
ip = 0...27T, for g = 1, where 7 RJ 0.308 is the triple point. For ip = we 
have 7 = 0.312 and the three states exist without complex continuation. The red 
square marks the ground state with real re. Its trace on the parameter space circle 
is represented by the red solid line. The two states with positive and negative 
imaginary parts of the eigenvalues re are represented by the blue triangle (blue 
dotted line for ip ^ 0) and the green circle (green dashed line), respectively. 





Figure 7. Exceptional point behaviour: circle with complex A = 0.04e lv , 
ip = . . . 27r, for g = l and 7 = 0.308. The triple point is located exactly at 
A = 0. The symbols and lines have the same meaning as in Fig. [5] The circle in 
the asymmetry parameter clearly reveals the nature of an EP3. 



3. The Matrix Model 

Here we present a simple low-dimensional matrix model that simulates the parameter 
dependence of the spectrum of the numerical solutions of the Gross-Pitaevskii equation 
as discussed in the previous section. In particular, the various types of the EPs as 
found above are of crucial importance: the simulation must reflect the particular types 
of these spectral singularities. 

3.1. General properties 

Our starting point is a two-mode model for a "PT-symmetric Bose-Einstein condensate 
[5] that mimics well that of the full Gross-Pitaevskii equation. In the two-mode 
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Figure 8. Eigenvalues E2, -B3 and E4 from Eq. 1(6} as a function of the non- 
Hermiticity parameter x for different strengths of the non-linearity, viz. c = 0.6 
(a), c = 0.1 (b). Note the similarity with the spectra in figure [4] 



approximation the 'P7 _ -symmetric stationary GP-equation assumes the form 
-i 1 + g\cb 1 \ 2 1 



The off-diagonal element in the Hamiltonian describes the tunnelling between the two 
modes, and has been normalised to unity. The analytic extension of the GP-equation 
in the spirit of Ref. [28] then yields the eigenvalues [6] 



Ex = - VI - x 2 , 



E 2 = + s/i 

= c - 



I 


— X 2 


-c 2 




X 2 + 


c 2 


I 


— X 2 


-c 2 




X 2 + 


c 2 



E A = c + x\ . (6) 

V ar + cr 

Here c = <?/2, and we use the symbol x instead of 7. By normalising the off-diagonal 
element to unity we have reduced the relevant interval of x between and 1. In 
figure [8] such spectra are illustrated for some values of c. It is striking how well 
these eigenvalues quantitatively describe the behaviour of those found numerically 
from the GP-equation. This perfect simulation includes in particular the spectral 
singularities and also allows to study the precise behaviour of the non-uniform limit 
j-)0 appearing here as c — > 0. 

At the critical point x = \[\ — c 2 the three eigenvalues E2 , E3 and E4 all assume 
the same value c. Beyond this point E3 and E4 become complex; in fact, E3 and £4 
have a square root singularity that appears to have some characteristics of a usual 
EP. The precise nature, especially its interplay with E2, cannot be deduced from the 
spectrum alone. It is the interaction giving rise to particular wave functions that 
determines the precise nature of a spectral singularity. As we are interested in the 
intersection point of the three levels E 2 , E 3 and E4 we restrict ourselves to a three- 
dimensional matrix model and leave out the level E\ altogether as it seems immaterial 
to our specific interest here. 

We seek a three-dimensional matrix which produces the spectrum {E2, £3, E4}, 
but in such a way that all levels are fully interacting. This is in contrast to the matrix 
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given in [B] where an interaction between the two branches E\, E 2 and E%, E4 is not 
taken into account. In other words, we seek a similarity transformation s such that 

h^s-j-s- 1 (7) 

where j is diagonal with the levels {E 2 , E3, E4} in its diagonal and s invokes the 
interaction between the levels. For our purpose the matrix s must bear the singularities 
of j. Moreover, a particular rank drop at the critical point x — \/l — c 2 determines 
as to whether the spectrum is simply degenerate, or whether it has a coalescence of 
two or of three levels, the former being a usual EP2 and the latter an EP3. Note that 
does not exist if s does not have full rank whereas h is still well defined but may 
no longer be diagonalisable; this is the signature of an EP. If rank(s) = l at the critical 
point, the Jordan normal form of h is expected to be 



J[h] = c 1 , (8) 




being the signature of an EP3: three levels are coalescing. 

A judicious choice fulfilling these requirements is given by 



1 y/i-x 2 v / r r 



2 



x'i+c> C+xy x . i+c .i _ (gj 



1 c- 



It can be checked that h yields Eq. (|8|) at x = V 1 — c 2 . Moreover, for x ^ vT— c 2 
one sees that encircling the critical point in the complex plane swaps the levels E3 
and E4, including their (unnormalised) eigenfunctions obtained from Eq. ([9]) 



/ VI - x 2 \ /VI 



'l_a;2_ c 2 

and I c ^ X 'V x 2 +c 2 



being just the characteristic of an EP2. In fact, the square root yl — x 2 — c 2 simply 
changes sign when a circle is traced out in the complex x-plane around the point 
VI c 2 thereby interchanging E3 and E4 and their corresponding eigenvectors. Note, 
however, that at the critical point all three eigenvectors also become equal (up to a 
constant) thus indicating the special nature of this point. 

We emphasise that the clear signature of an EP3 as presented above is only 
available when analytic expressions are at hand. To obtain this result directly by 
numerical means is virtually impossible. One rather has to resort to the means as 
presented in the previous section, where the existence of the EP3 is established by 
invoking an asymmetry and then use closed contours in a suitable complex parameter 
plane. To elucidate this pattern we discuss in the appendix a three-dimensional 
matrix demonstrating the underlying mechanism; as our matrix h is unsuitable for 
demonstration due to its complicated structure we use a much simplified matrix. 

3.2. Behaviour for c — > 

All these statements hold for c > but care must be taken when c — > 0. Recall that 
the GP-equation has a non-uniform limit for g — ¥ as is discussed in Section 2. The 
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same holds for the matrix model when c — > 0. As the limiting behaviour is subtle 
and difficult to obtain numerically, we here give some results obtained in the matrix 
model. Note in particular how the triple point approaches the branch point where the 
■pT-eigenfunctions meet when c — > for the matrix model (figure [5} and g — > for 
the numerical solution (figure SJ. 

The eigenvalues E3 and E4 have, for c > 0, the value c at x = 0. Yet these 
eigenvalues are ±1 for c = and x = 0. Similarly, the derivative tends to infinity at 
x — for c — > 0, as can be checked analytically, and is visualized also in figured] But 
it is zero when c = is taken from the outset. In other words, E4 is getting nearer 
and nearer to E%, but at x = in a non-continuous fashion (E3 approaches — E4 
for c — > 0). Note that this behaviour nicely reflects the behaviour of the numerical 
solutions of the full GP-equation (see figure 0]). 

At x = 1 it is the eigenfunctions that show the typical non-uniform behaviour. 
Note that the eigenvectors as listed in Eq. © cannot be used in the limits x — > 1 
and c — > 0, irrespective of the order of the limits taken, since s _1 does not exist. 
Rather, the full matrix h must be considered. If c = 0, corresponding to switching 
off the non-linearity in the GP-equation, the eigenvectors of the two positive solutions 
(E2 = E4 = \/l — x 2 ) are given by 



whereas, when the limits are commuted (x — 1 first and then c — > 0) the corresponding 
eigenvectors are 



Recall that E3 and E4 are complex when x is to the right of the critical point, that 
is if x > \/l — c 2 , a condition met for the order of the limits taken. This is reflected 
here in the complex eigenvectors, even though the eigenvalues are zero in either limit. 

A further point of interest is the behaviour of the scalar product taken from the 
eigenvectors referring to the plain solution for c = and the eigenvector referring 
to c > 0. Of interest is the limit c — > of the scalar product. For the eigenvalue 

E4 = c + xJ 1 ~i^_~£ 2 we obtain for the normalised eigenvectors the scalar product 



We simply list some major facts as they can be easily verified from Eq. (|10[) . The 
scalar product behaves in its limit c — > non-uniformly at both ends of the a;-interval. 
At x = Eq. (fit))) assumes the value (1 + c) + Ac 2 whereas the value -y/2/3 is 
assumed when c = is set first. At the right end of the x-interval the value is unity 
when c = while the value is \J (1/2 + i) (3 — i) /5 when x = 1 is set before c — > 0. 
We mention that the intricacy of the expansion in powers of c probably hints to the 
difficulty in approaching the limit c — > numerically using the GP-equation directly. 
In fact the expansion reads 






(10) 



2V2> 2 - 1) 



+ 



V3 - 4a; 2 -I- x 4 



{x 2 - 3)V3-4x 2 + .t 4 
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Figure 9. Scalar product of the eigenvector for c = with the eigenvector 
for c' ^ as function of x in the range Q < x < 1 (a). The values displayed 
refer to c' = 0.1 (red solid line), c' = 0.05 (green dashed line), c' = 0.01 (blue 
dotted line) and c' — > (purple dashed-dotted line). For better illustration the 
x-axis is stretched at x = (b) and at X = 1 (c). All lines are purely real only 
for x < Vl — c 2 ; it explains the ending for x < 1 as only pure real values are 
displayed. 



where the x- dependent coefficient of yfc vanishes for x — 1; it means that the y/c- 
behaviour will show for x < 1 but not for x = 1. In figure [9] the scalar product 
is illustrated for a few values of c. The non-uniform behaviour at both ends of the 
x-interval is clearly visible. We note that the scalar product assumes the value y / 2/3 
at the critical point, independent of c. 

4. Summary and Conclusions 

A number of novel aspects originating from the non-linear term g\ip\ 2 in the interaction 
term of the Schrodinger equation are communicated in the present paper. There is the 
occurrence of two new states in the low energy spectrum. Naively, they are seen for 
g > only when they emerge as complex solutions from the ground state for 7 > 7 cr 
thereby breaking the "PT-symmetry of the system. An appropriate complex extension 
of the GP-equation reveals that two real additional solutions do exist for 7 < 7 cr 
which are in fact eigenfunctions of VT- A further novel result is the triple point at 
7 = 7 cr , in particular the singular spectral behaviour. With a superficial look the 
singularity appears as a "common" EP2 where the two new solutions turn complex 
while the ground state from which they emerge seems to be unaffected. However, a 
more thorough analysis reveals that the three states are nested in an EP3; in fact, for 
any generic perturbation this result prevails. 

While these results are obtained by careful numerical analysis the point g = 
being a singular point of the non-linear equation seems to be inaccessible by numerical 
means. In other words, the limit g — > leading to the simpler linear equation is subtle, 
there is a clear non-uniform behaviour for the spectrum and eigenfunctions. To study 
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this aspect in analytic terms a matrix model simulating quantitatively the spectrum 
of the GP-equation is presented. While the model elucidates the limit very clearly it 
also sheds light on the specific behaviour of the triple point: the distinction between 
the EP2 and EP3 behaviour depending on the specific approach of the singular point 
becomes evident in analytic terms. The matrix model is designed to be as simple 
and instructive as possible since our interest is focussed upon the explanation of the 
spectra. An exact quantitative agreement with the numerical results is not intended 
as it will complicate the analytic discussion. 

We believe that these findings are expected to have an impact in experimental 
work in that specific effects can be detected related to the structure of the model 
studied. It is well known that EPs have a definite physical meaning [31 in a 
great variety of physical systems. The presence of the triple point has already been 
established to lead to an instability in the dynamical behaviour [S]. It remains a 
challenge to detect the signature of the EP3 in a BEC. Whether or not even the two 
additional states have physical significance is, at this stage, an open question. 
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Appendix 

The matrix 




2 -1 
c + x -1 
-1 

has the spectrum 

Ei = 1 

E-2 = 1 + \J\-c-x 
E 3 = l- Vl-c-x. 

Clearly, and from the structure of the associated eigenfunctions, there is an EP2 at 
x = 1 — c. However, when x = 1 — c is chosen from the outset the Jordan normal form 
reads 

110' 
1 1 
1 

which clearly indicates an EP3. How to find this by numerical means? 
We use a perturbation e 

2 + e -1 
2 + c + x —1 — c — x 
-10 2 
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The expansions of the eigenvalues read 

£i = l + - e + 0(e 2 ) 

1 — c — x 

E 2 = l + ^l-c-x+ ( \ +C + \ + n 7= 1 )e + Q(e 2 ) 
\2(l-c-x) 2vT -c-xj 

E 3 = l- VT^c—x + ( \ +C + X . 0/1 1 ) e + 0( £ 2 ). 

\2(l-c-a;) 2vl-c-a;/ 

Obviously, the expansion breaks down at the critical point x = 1 — c. When x is set 

equal to 1 — c in the perturbed matrix, we now see clearly how the eigenvalues sprout 

from the EP3 in the expected manner. The expansions read now 

£?i = 1 - 2 1/3 cxp(2i7r/3)e 1/3 + 0(e 2/3 ) 

E 2 = l-2 1 /V/ 3 + 0(e 2 / 3 ) 

£ 3 = 1- 2 1 / 3 cxp(-2iV3)e 1/3 + 0(e 2 / 3 ). 
This pattern is basically unchanged for any generic perturbation. It clearly 
demonstrates the principle underlying in the numerical manifestation of the EP3. 
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